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1 . INTRODUCTION 

As  the  development  of  sonar  systems  becomes  increasingly  sophisticated,  so 
the  need  for  reliable  undersea  acoustic  modelling  becomes  more  acute.  Particul- 
arly in  the  case  of  low  frequency  long  range  propagation,  environmental  para- 
meters such  as  sound  speed  profiles  and  bathymetry  are  likely  to  exhibit  signifi- 
cant changes  over  the  range  of  interest.  For  such  cases  predictions  based  upon 
range-independent  models  become  dubious,  and  an  effective  range-dependent  model 
is  essential. 

There  are  a number  of  approaches  to  the  problem  of  range  dependence  including 
variations  and  combinations  of  ray-tracing  and  normal  mode  techniques  (see  for 
example,  references  1,2,  and  3),  but  it  is  still  true  to  say  that  no  single  model 
is  ideal  for  all  cases.  One  of  the  most  significant  milestones  in  the  recent 
history  of  range-dependent  modelling  has  been  the  application  of  the  parabolic 
equation  technique  to  underwater  acoustics(ref .4) . This  has  given  rise  to  what 
could  be  described  as  the  first  range-dependent  direct  field  propagation  model. 

This  report  describes  a range-dependent  acoustic  propagation  loss  model  based 
upon  the  parabolic  equation  technique  and  developed  for  use  on  the  DRCS  IBM  370/168 
computer.  An  earlier  version  of  the  model,  derived  from  one  in  use  at  SACLANT 
ASW  Research  Centre(ref .5)  has  been  implemented  by  the  author  on  a UNIVAC  1108 
computer  at  U.S.  Naval  Underwater  Systems  Center,  New  London  (NUSC ) , (ref .6) . 

The  new  version  described  in  this  report  differs  greatly  from  the  version  in 
use  at  NUSC.  The  starting  field  is  computed  directly  rather  than  by  invoking 
an  elaborate  and  often  extremely  time  consuming  normal  mode  calculation  procedure. 
The  second  main  difference  is  in  the  handling  of  range  variation  of  parameters, 
and  the  treatment  of  bottom  modelling.  The  new  model  assumes  continuous  inter- 
polation in  range  between  successive  sound  speed  profiles,  giving  a greatly 
improved  capability  for  handling  such  phenomena  as  the  reflection  of  energy  from 
sloping  bottoms.  Like  other  models  in  common  use,  the  parabolic  equation  model 
works  extremely  well  within  the  limits  for  which  it  is  intended,  but  cannot  be 
expected  to  yield  reliable  results  if  those  limits  are  exceeded.  The  point  here 
is  that  the  model  will  still  produce  answers,  and  the  answers  will  still  look 
reasonable,  but  they  will  nevertheless  be  systematically  incorrect.  This  is 
potentially  dangerous  as  an  inexperienced  user  not  following  the  guidelines  may 
not  recognise  when  his  predictions  are  useless.  One  of  the  purposes  of  this 
report  is  to  explain  the  mechanics  of  operation  of  the  model.  At  the  same  time 
it  is  intended  that  the  report  should  provide  an  understanding  of  the  reasoning 
behind  the  various  procedures  in  the  model  so  that  a user  may  better  recognise 
when  the  model  is  or  is  not  suited  to  a particular  application,  and  may  use  the 
model  with  maximum  efficiency. 

An  appendix  describing  the  input  data  format  is  included  in  this  report. 

Copies  of  the  program  listing  are  available  from  the  author  on  request. 


2.  THE  PARABOLIC  EQUATION  AND  THE  MARCHING  SOLUTION 

The  following  section  merely  outlines  the  main  steps  leading  to  the  marching 
solution  of  the  parabolic  approximation  to  the  time  invariant  reduced  wa.e 
equation  (Helmholtz  equation).  For  a more  complete  derivation  the  reader 
should  consult  references  4,5  or  6. 

The  Helmholtz  equation  may  be  approximated  by  a parabolic  equation  of  the  form 
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where 


ip  = is  a wave  function 

kQ  = reference  wave  number  (=  w/c0) 

r i = index  of  refraction  (=  c /c) 

o 


c = sound  speed  c = c(r,z) 
a)  = source  angular  frequency 

The  main  assumptions  are  those  of  cylindrical  symmetry  and  the  paraxial 
(small  angle)  approximation 


2k0  9 £ 
8r 


(2) 


which  is  satisfied  by  restricting  the  wave  normals  in  the  field  to  small  grazing 
angles.  Experience  shows  that  the  model  gives  quite  good  results  when  the  wave 
normal  grazing  angles  are  restricted  to  less  than  about  ± 25°.  Snell's  law  may 
be  used  to  determine  the  range  of  source  grazing  angles  corresponding  to  25°  wave 
normal  grazing  angles  at  sound  speed  profile  minimum. 

Solution  of  equation  (l ) involves  two  additional  assumptions.  The  first  is 
that  both  ip  and  dp_  go  smoothly  to  zero  for  z ’*-<*>.  This  is  achieved  by 
dz 

extending  the  field  far  enough  into  the  bottom  in  the  presence  of  bottom  atten- 
uation, and  reflecting  the  extended  field  about  the  sea  surface  with  the  approp- 
riate 180°  phase  shift.  The  second  assumption  is  that  at  any  range  the  refrac- 
tive index  U is  independent  of  range  and  depth.  While  this  is  obviously  incorrect 
in  any  range  dependent  problem  of  interest,  the  errors  introduced  by  the  assump- 
tion may  by  small  enough  for  useful  results  to  be  obtained.  It  can  be  shown 
that  this  assumption  paves  the  way  for  a marching  solution  of  equation  (l ) of  the 
form 

<fi(r  + A r,z)  = exp  (77 2 _ 1 ) Ar).  3 { exp  (~^r  S2)  3-  \ip{r,z)\  1 (3) 

where  3 and  3 ^ denote  the  forward  and  inverse  Fourier  transforms,  respectively. 

Beginning  with  a known  field  at  one  range  ^ (rQ,z),  the  recursive  form  of 

equation  (3)  enables  us  to  find  the  field  at  all  other  ranges  by  "marching"  the 
solution  out  in  range.  Use  of  a Fast  Fourier  Transform  (FFT)  algorithm  enables 
this  to  be  carried  out  extremely  rapidly. 


3.  THE  STARTING  FIELD 

One  of  the  big  problems  with  generation  of  a starting  field  by  means  of 
computation  of  normal  modes  has  been  that  in  the  large  proportion  of  cases  of 
interest,  calculation  of  some  of  the  contributing  modes  is  difficult,  and  up  to 
several  hundred  modes  may  be  required.  Double  or  extended  precision  arithmetic 
is  usually  necessary  and  computing  times  are  often  lengthy.  While  such  an 
approach  has  been  fairly  standard  for  initiating  parabolic  equation  models  up 
until  now  a careful  consideration  of  the  problem  shows  that  the  precision  neces- 
sary for  determination  of  the  modes  is  far  greater  than  that  necessary  to 
describe  the  starting  field  itself,  and  that  the  normal  mode  method  is  grossly 
inefficient. 
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A number  of  laboratories  using  parabolic  equation  models  have  realised  this 
problem,  and  currently  a number  of  methods  are  being  investigated  to  overcome 
it (ref .7). 

The  starting  field  generator  used  for  the  presently  described  model  is  an 
extremely  simple  one,  made  possible  by  the  existing  constraints  on  the  model  of 
dependence  of  FFT  size  on  frequency  and  the  paraxial  limitation.  This  generator 
consists  of  a simple  algorithm  for  calculating  a table  of  phase  and  amplitude 
values  at  a series  of  equispaced  points  across  a vertical  line  a short  distance 
from  the  location  of  a point  source.  For  the  purposes  of  generating  the  initial 
field,  sound  speed  is  assumed  to  be  locally  depth  - invariant  at  the  source,  so 
that  the  wavefronts  intersecting  the  line  may  be  considered  to  be  spherical. 

The  table  values  must  be  calculated  at  suitably  close  depth  intervals  along  the 

line  for  a distance  subtending  the  source  angle  of  interest,  with  table  values 

outside  this  angle  set  to  zero.  The  paraxial  approximation  limitation  of  the 

model  precludes  source  grazing  angles  much  in  excess  of  25°  which  in  turn  restricts 

the  vertical  coverage  of  the  non-zero  data  values  (referred  to  as  the  aperture) 

and  so  limits  the  error  introduced  by  considering  that  the  source  region  sound  > 

speed  is  depth  invariant. 

Using  the  analogy  that  the  aperture  behaves  like  a beamforming  array,  it  is 
apparent  that  for  satisfactory  performance  of  this  starting  field  model,  two 
conditions  must  be  satisfied.  The  first  condition  is  that  the  spatial  sampling 
be  no  coarser  than  half  wavelength  to  avoid  aliasing  problems.  The  second  is 
that  the  aperture  should  be  sufficiently  large  that  it  does  not  impose  significant 
diffraction  effects  upon  its  model  of  the  source.  We  will  deal  with  the  spatial 
sampling  first.  The  comments  apply  equally  well  to  the  rest  of  the  field  as  to 
the  starting  field. 

3.1  Choice  of  FFT  size 

The  spatial  sampling  is  finer  than  half-wavelength  when  the  minimum  size 
of  the  FFT  (the  number  of  complex  data  points  NFFT)  appropriate  to  the 
problem  is  simply  the  smallest  integer  power  of  two  satisfying  the  inequality 

4D  A.  min  /.  \ 

NFFT  2 u ' 

which  may  be  arranged  to  a more  useful  form 

NFFT  ? ; 8D  (2) 

Amin  ' ' 


where  D is  water  depth  and  X min  is  the  wavelength  associated  with  the  sound 
frequency  under  consideration  at  the  sound  speed  profile  minimum.  The  factor 
of  four  in  (l  ) arises  from  the  method  of  satisfying  the  boundary  conditions 
noted  in  Section  2.  The  field  is  extended  into  the  bottom  a distance  equal 
to  the  water  depth,  and  this  extended  field  is  then  reflected  in  the  sea 
surface,  phase  reversed,  requiring  that  the  model  handle  a spatial  field  four 
times  greater  than  that  necessary  to  describe  the  water  column  alone. 

To  give  an  example,  consider  100  Hz  propagation  in  water  5000  m deep. 

If  the  minimum  sound  speed  is  1450  m/s“  , A min  will  be  14.5  m at  100  Hz,  so 
that  from  (2)  we  require  NFFT  £ 2759.  A 4096  point  trahsform  is  the  small- 
est satisfying  this  requirement.  Experience  with  use  of  the  model  however 
shows  that  for  cases  in  which  only  low  grazing  angle  energy  is  of  interest, 
it  would  be  possible  to  economize  on  computing  by  using  only  a 2048  point 
transform  sampling  slightly  in  excess  of  A min. 

2 

The  explanation  lies  in  the  application  of  sampling  theory  to  digital 
transforms  such  as  the  FFT.  The  Fourier  transform  for  a spatially  sampled 
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waveform  over  a vertical  aperture  is  the  plane  wave  angular  spectrum  which 
would  result  from  either  an  infinitely  long  aperture  or  a continuum  of  sample 
points  over  the  aperture,  convolved  with  a Dirac  comb  function  where  the 
angular  spacing  of  the  comb  teeth  is  given  by  (e.g.  reference  8) 

Sin  0 = An  / \ 

z 


where  ^ is  the  (constant)  wavelength  in  a medium  of  constant  sound  speed, 

N is  the  total  number  of  evenly  spaced  samples  of  the  waveform  over  the 
vertical  aperture,  and  z is  the  vertical  extent  of  the  aperture.  In  the 
present  context  N is  equivalent  to  NFFT,  z is  4D,  and  the  constancy  of  A 
is  one  of  the  basic  approximations  of  the  parabolic  equation  technique. 

Since  the  right  hand  side  of  (3)  will  often  be  greater  than  one,  it  is 
apparent  that  the  angular  spacing  of  the  comb  function  may  be  complex. 

Indeed  the  fundamental  spectrum  itself  may  include  complex  wavenumbers , but 
in  the  present  application  most  of  the  energy  lies  in  the  real  portion  of  the 
spectrum  and  the  complex  portion  may  be  ignored.  Provided  that  AN?  2 

z 

there  will  be  no  overlap  of  the  fundamental  angular  spectrum  with  the  complex 
higher  order  angular  spectra,  so  that  they  may  for  practical  purposes  be 
ignored.  As  An  is  allowed  to  fall  below  2 however  (in  other  words,  the 
z 

sample  spacing  is  permitted  to  exceed  half  wavelength),  then  some  overlap  of 
the  various  spectra  will  occur.  Where  the  half  wavelength  spacing  is  ex- 
ceeded only  slightly,  the  overlap  also  will  be  slight,  with  large  positive 
angles  of  the  fundamental  spectrum  overlain  by  that  portion  of  the  adjoining 
higher  order  spectrum  representing  large  negative  angles,  and  vice  versa. 

This  phenomenon  may  be  referred  to  as  angular  aliasing. 

The  reason  that  we  can  tolerate  some  apparent  undersampling  of  the  vertical 
field  distribution  now  becomes  clear.  Since  the  parabolic  equation  model  is 
limited  by  considerations  of  accuracy  to  * 25°  grazing  angles,  then  the  energy 
in  the  source  function  angular  spectrum  should  be  essentially  confined  within 
these  limits.  If  higher  angles  in  the  source  spectrum  have  no  energy  assoc- 
iated with  them,  then  provided  that  energy-free  aliased  angles  merely  overly 
other  energy-free  high  angles,  there  will  be  no  degradation  of  the  - 25° 
portion  of  the  spectrum  due  to  the  apparently  coarse  spatial  sampling. 

To  take  full  advantage  of  this  it  is  clear  that  we  can  relax  the  constraint 
on  sampling  of  the  inequality  (2).  For  a limiting  source  angle  y , an  alias- 
ing folding  angle  exists  such  that  when  the  spectrum  is  folded  about  it,  the 
aliased  portion  of  the  spectrum  will  not  overlap  the  source  angle  limits. 

This  condition  is  satisfied  when  (modifying  (2)),  we  have 

NFFT  £ 4 (l  + Sin  v)  (4) 

Vmin 

Returning  to  our  example,  if  y = 25°,  we  find  from  (4)  that  NFFT  £1955, 
rather  than  2759  as  implied  by  (2),  so  that  a 2048  point  transfrom  is  indeed 
adequate  in  this  instance. 

In  the  present  version  of  the  model  the  largest  available  transform  is  one 
of  4096  points,  which  is  satisfactory  for  most  cases  of  interest.  This 
could  be  simply  extended  by  increasing  array  dimensions  provided  that  a 
suitable  FFT  algorithm  is  used.  The  MKLFFT  routine  originally  used(ref.9) 
is  limited  to  a maximum  of  4096  points  without  modification,  but  a faster 
machine  language  FFT  routine  now  in  use  at  DRCS  can  handle  arrays  up  to  2^9. 
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3.2  Starting  field  aperture 

Choice  of  aperture  size  for  the  starting  field  generator  is  a compromise 
between  opposing  factors.  The  model  for  simplicity  assumes  a constant 
sound  speed  over  the  starting  field  aperture.  This  will  generally  be  a 
good  approximation  over  a small  aperture,  and  a poor  one  for  a large  aperture. 
The  aperture  should  thus  be  made  no  larger  than  necessary  to  minimise  potent- 
ial errors  arising  from  this  assumption.  On  the  other  hand,  the  aperture 
should  be  large  enough  that  diffraction  effects  associated  with  it  are  small 
compared  with  its  geometrical  properties  if  it  is  to  effectively  represent 
the  source. 

Since  by  virtue  of  the  application  of  Fourier  transform  techniques  to  the 
problem  we  are  resolving  the  field  into  its  plane  wave  components,  we  may 
therefore  resort  to  simple  diffraction  theory  for  analysing  the  behaviour 
of  the  aperture.  Figure  1 shows  the  variation  of  the  half-angle  between 
3 dB  points  for  the  main  beam  of  a Fraunhofer  diffraction  pattern  expressed 
as  a function  of  the  aperture  length  in  wavelengths.  Substantial  reductions 
in  beamwidth  are  effected  by  modest  increases  in  aperture  length  out  to 
about  1 0 A , and  certainly  the  aperture  should  if  possible  be  no  smaller  than 
5*.  For  a limiting  source  angle  of  15°  or  larger,  an  aperture  of  5 A.  or 

more  should  not  seriously  distort  the  starting  field.  However,  if  we  wished 

to  consider  a limiting  source  angle  of  only  5 , then  it  is  apparent  from 
figure  1 that  the  minimum  aperture  length  should  be  increased  to  nearer  20.  . 

For  programming  simplicity  we  will  arbitrarily  choose  an  aperture  10  A in 
length.  The  source  to  aperture  distance  will  then  be  determined  simply  by 
the  limiting  source  angle,  which  under  the  circumstances  may  lie  roughly 
between  10°  and  25°. 

For  the  25°  limiting  source  angle  case,  the  source  to  aperture  distance 
will  be  about  10A  also.  The  assumption  of  constant  sound  speed  will  intro- 
duce phase  errors  into  the  starting  field,  but  for  normally  encountered 
sound  speed  gradients  these  phase  errors  will  be  quite  small.  For  an  in- 
water sound  speed  gradient  of  0.2  s-^ , and  a frequency  of  100  Hz,  the 

assumption  of  constant  sound  speed  will  amount  to  a maximum  error  in  sound 
speed  over  a 10A  aperture  of  the  order  of  This  in  turn  will  yield  a 

maximum  phase  error  over  the  source  to  aperture  distance  equivalent  to  about 
A/10,  which  corresponds  to  an  uncertainty  in  the  position  of  the  source  of 
only  a few  metres.  While  the  phase  errors  due  to  the  assumption  of  locally 
constant  sound  speed  are  proportional  to  frequency,  the  aperture  is  scaled 
to  wavelength,  which  tends  to  compensate  for  this.  For  cases  where  large 
apertures  are  necessary  because  of  very  low  frequencies  the  average  sound 
speed  gradient  across  the  aperture  will  typically  be  considerably  less  than 
0.2  s-',  so  that  for  most  practical  cases  the  errors  introduced  by  this 
simple  source  model  will  be  quite  insignificant. 

It  should  be  emphasized  that  the  approximation  of  locally  constant  sound 
speed  over  the  aperture  is  only  invoked  for  the  calculation  of  the  starting 
field.  Once  the  starting  field  has  been  determined  the  field  is  then 
calculated  at  all  successive  range  intervals  by  the  marching  technique  of 
equation  (3),  using  the  full  sound  speed  profile.  The  predictions  thus 
generated  with  the  simple  starting  field  technique,  provided  that  the  above 
constraints  on  aperture  size  are  observed,  will  be  virtually  indistinguish- 
able from  those  obtained  using  an  equivalent  normal  mode  starting  field 
generator. 


4.  SOUND  SPEED  RANGE  VARIATION  AND  BOTTOM  MODELLING 

Sound  speed  profiles  consisting  of  up  to  100  depth,  speed  pairs  must  be  read 
in  for  both  zero  and  maximum  ranges  from  the  source  as  well  as  at  any  desired 
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intermediate  ranges.  The  bottom  profile  point  is  interpreted  as  being  at  the 
sea  bottom.  A maximum  in-bottom  sound  speed  is  read  in,  as  well  as  a chosen 
in-bottom  sound  speed  gradient.  Both  of  these  parameters  apply  over  the  full 
range  of  the  run. 

Each  input  profile  is  first  linearly  interpolated  in  depth  to  give  a profile 
of  points  at  the  spatial  frequency  of  the  FFT.  At  every  range  step  a new  local 
sound  speed  profile  is  horizontally  linearly  interpolated  between  the  two  applic- 
able interpolated  input  profiles.  The  bottom  depth  for  the  intermediate  profile 
is  also  linearly  interpolated.  Where  one  of  the  input  profiles  has  reached 
bottom,  in  the  case  of  a sloping  bottom,  the  horizontal  interpolation  procedure 
must  be  modified.  If  the  water  depth  is  greater  than  1000  m,  the  sound  speed  is 
chosen  to  be  the  same  as  the  sound  speed  at  that  depth  f or  the  profile  which  has 
not  yet  reached  bottom.  Where  the  water  is  less  than  1000  m deep,  the  inter- 
mediate sound  speed  is  interpolated  between  the  in-water  value  at  the  appropriate 
depth  on  the  deeper  side,  and  the  in-water  bottom  value  on  the  shallower  side. 

The  break-point  value  of  1000  m is  chosen  as  the  nominal  depth  of  the  bottom 
of  the  thermocline,  below  which  the  range  dependence  of  sound  speed  profile  can 
be  expected  to  be  very  slight  in  comparison  with  the  variations  closer  to  the  sea 
surface. 

At  every  range  step  the  profile  is  extended  into  the  bottom  by  increasing  the 
in-water  bottom  sound  speed  value  according  to  the  chosen  input  sub-bottom  sound 
speed  gradient  until  the  maximum  sub-bottom  sound  speed  is  reached,  whereafter 
the  speed  will  be  constant. 

This  method  of  profile  interpolation  with  range  has  the  effect  of  smoothing 
out  the  in-water  range  variations  into  a large  number  of  small  increments,  as 
distinct  from  having  relatively  large  discontinuities  in  profile  at  a few  ranges. 
The  smoother  approach  must  give  a more  realistic  description  of  what  occurs  in 
the  sea  in  most  situations,  and  certainly  is  better  suited  to  the  capabilities 
of  the  model(ref .6) . 

4.1  Choice  of  bottom  gradient 

The  interaction  of  acoustic  energy  with  the  sea  bottom  is  a source  of 
difficulty  in  all  acoustic  models,  but  especially  those  specifically  suited 
to  dealing  with  low  frequencies.  Before  looking  at  the  present  model  we 
will  consider  two  other  approaches. 

In  a raytracing  model  it  is  usual  to  determine  the  proportion  of  energy 
reflected  from  the  bottom  by  referring  to  tables  of  bottom  loss  versus 
grazing  angle  for  various  frequencies.  Such  tables  are  generally  determined 
experimentally,  and  different  tables  may  apply  to  different  regions. 

Figure  2 shows  typical  experimental  data. 

The  Fast  Field  wave  model ( ref .1 0)  does  not  use  a bottom-loss  table. 

Instead  the  bottom  is  modelled  as  comprising  one  or  more  discrete  fluid  layers, 
and  the  layer  sound  speeds  and  densities  are  varied  to  produce  various 
combinations  of  plane'  wave  reflection  coefficients  as  a function  of  grazing 
angle.  Figure  3 shows  some  bottom  loss  predictions  calculated  on  the  basis 
of  a simple  two  layer  plane  wave  reflection  coefficient  model,  for  likely 
combinations  of  bottom  and  sub-bottom  sound  speeds  and  densities.  Where 
the  lateral  wave  is  taken  into  account(ref . 1 1 ) or  a multilayer  bottom  is 
assumed,  such  smooth  curves  will  generally  be  replaced  by  a more  complicated 
interference  pattern  set  up  between  multiple  bottom  arrivals. 

Multilayer  bottom  modelling  is  generally  only  significant  at  low  frequen- 
cies where  sub-bottom  attenuation  is  low  enough  for  significant  energy  to  be 
returned  to  the  water  from  deep  layers  and  where  sufficient  information 
about  the  sub-bottom  structure  exists  to  give  the  model  some  practical  basis. 
Figure  4,  following  reference  10,  is  an  example  of  a multilayer  bottom  loss 
prediction. 
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The  parabolic  equation  model  uses  neither  of  the  above  approaches. 

Since  it  is  a field  method  it  does  not  keep  track  of  individual  rays,  so  the 
first  method  is  inapplicable.  The  second  method  is  not  suitable  either, 
since  the  parabolic  equation  model  is  inherently  incapable  of  accommodating 
a discontinuity  in  the  sound  speed  profile. 

Instead  we  model  the  bottom  interface  by  a combination  of  a relatively 
sharp  gradient  in  sound  speed  and  a volume  absorption  factor.  The  gradient 
causes  some  energy  to  be  refracted  back  into  the  water  before  it  has  suffer- 
ed too  much  absorption.  Control  over  the  amount  of  energy  refracted  up- 
wards is  established  by  varying  the  sub-bottom  gradient  and  the  absorption. 
To  understand  the  way  in  which  the  parabolic  equation  bottom  loss  modelling 
works  it  is  instructive  to  picture  what  happens  in  terms  of  ray  theory. 

Consider  energy  impinging  on  the  bottom  to  be  represented  by  a ray  equiv- 
alent, incident  on  the  bottom  at  an  angle  to  the  normal  6q,  with  a sound 

speed  cq.  From  Snell's  Law 


p = Sin  0 /c 
o o 

where  p is  known  as  the  ray  constant.  If  we  represent  the  bottom  by  a 
region  with  constant  sound  speed  gradient 


dc 

dz 


g 


we  have  (following  reference  12) 


where  <P  is  the  local  ray  angle,  S is  the  distance  along  the  ray,  and  p is 
the  radius  of  curvature. 

This  simply  states  that  the  curvature  of  the  ray  path  in  the  presence  of 
a constant  gradient  is  constant,  that  is  to  say,  the  path  is  circular. 

It  is  apparent  from  figure  5 that  the  ray  will  emerge  back  into  the  water 
from  the  sub-bottom  region  of  constant  gradient  at  the  same  angle  as  it 
entered,  reflected  about  the  normal.  The  path  length  of  the  ray  in  the 
sub-bottom  region  is  given  by 

s = ( 7i  - 2 0Q)p 


Figure  6 shows  the  variation  of  sub-bottom  distance  S(gQ,  c ) travelled 

by  the  ray  equivalent  before  re-emergence  into  the  water  column,  as  a 
function  of  grazing  angle  for  the  case  of  a bottom  interface  soijind  speed 
cq  = 1500  m.s-^  and  a sub-bottom  sound  speed  gradient  gQ  = 1 s 

For  other  cases  the  sub-bottom  distance  may  be  simply  scaled  from  figure 
6,  using  the  relation 

S(g,  c)  = S(g  , c ) x c/ ( 1 500  g) 
o o 


(1) 
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where  c is  the  bottom  interface  sound  speed  and  g is  the  sub-bottom  sound 
speed  gradient. 

For  most  situations  of  interest  the  variation  of  c will  only  be  small, 
and  the  main  parameter  of  concern  will  be  the  sub-bottom  sound  speed  gradient 
g.  The  higher  the  value  of  g,  the  greater  will  be  the  curvature  of  the  sub- 
bottom ray  equivalent.  This  in  turn  will  result  in  a smaller  total  phase 
change  within  the  bottom.  Hard  bottoms  will  thus  be  represented  by  high 
values  of  g,  and  softer  bottoms  by  lower  values  of  g. 

A value  of  1 s is  a reasonable  choice  of  g for  a soft  bottom.  For  a 
hard  bottom  the  highest  value  of  g possible  is  obviously  going  to  provide 
the  best  approximation  to  a discontinuous  step  in  sound  speed.  However, 
the  higher  the  value  of  g,  the  shorter  is  the  permissible  range  step  for  the 
marching  solution,  and  consequently  the  greater  the  amount  of  computation 
needed  to  predict  propagation  loss  to  any  given  range. 

For  low  frequency  applications  there  seems  to  be  little  point  in  increas- 
ing g beyond  about  10  s-  , except,  perhaps,  for  shallow  water,  hard  bottom 
modelling  as  in  Section  6.  A comparison  of  a c.w.  prediction  for  bottom 
limited  propagation(ref .6)  using  the  FFP  model  with  a parabolic  equation 
prediction  assuming  a g of  4 s-^  shows  almost  identical  behaviour  of  the 
interference  fringes  arising  from  bottom  reflections  (figure  7).  One  could 
not  expect  the  two  curves  to  be  indistinguishable,  since  apart  from  the  in- 
herent approximations  of  the  parabolic  equation  technique,  the  two  methods 
treat  the  bottom  differently  and  interpolate  the  input  sound  speed  profile 
points  differently.  Nevertheless,  the  general  features  of  each  curve 
correspond  remarkably  well.  From  0-8  kyds  each  curve  shows  a rapid  shadow- 
zone  fall  off,  from  10-20  kyds,  a coarse  interference  structure  (which  corres- 
ponds to  singly  bottom-reflected  Lloyd's  mirror  fringes),  and  from  20-25  kyds 
a zone  of  fine  fringes  which  become  coarser  from  25-33  kyds. 

Within  the  full  spectrum  of  acoustic  propagation  prediction  models  only  the 
FFP  and  normal  mode  models  could  be  expected  to  produce  results  more  similar. 
When  the  parabolic  equation  model  with  its  additional  capability  for  handling 
range  dependent  data  also  shows  close  agreement  in  a range-independent  case 
with  such  a mathematically  powerful  model  as  the  FFP,  it  is  clear  that  the 
parabolic  equation  method  too  is  a powerful  technique. 

It  has  proved  far  easier  in  a study  of  bottom  limited  propagation  to  achieve 
close  correspondence  between  predictions  from  a number  of  models  than  between 
any  of  the  models  and  experimental  data,  despite  exceptionally  precise  experi- 
mental control(ref . 1 3) • At  present  it  is  still  usual  to  describe  the 
bottom  interaction  rather  like  an  electronic  "black  box",  as  though  all  effects 
take  place  at  the  boundary.  This  is  demonstrably  an  oversimplified  approach 
for  low  frequency  propagation  modelling,  but  is  probably  adequate  for  most 
cases  at  present  because  of  other  limiations,  including  navigational  error  and 
lack  of  precise  environmental  data. 

4.2  Sub-bottom  volume  absorption  coefficient 

It  should  be  noted  that  the  expression  S(g,c)  in  (l)  is  only  the  sub- 
bottom  distance  of  the  ray  equivalent.  To  obtain  the  bottom  loss  a the 
appropriate  angle  it  is  necessary  to  multiply  this  by  the  bottom  attenuation 
coefficient  a (dB/m). 

In  the  program  the  chosen  value  of  a is  simply  used  to  determine  the  mag- 
nitude of  an  imaginary  component  of  the  square  of  the  refractive  index  which 
appears  in  equation  (3).  This  then  causes  appropriate  attenuation  of  the 
sub-bottom  component  of  the  field. 

In-water  volume  absorption  could  be  treated  in  a similar  manner,  but  is 
included  in  the  present  model  instead  as  a simple  frequency  dependent  atten- 
uation coefficient  acting  upon  range  and  based  upon  the  empirical  formula  of 
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Thorp(ref . 1 4) . 

If  experimental  data  is  available  for  the  area,  the  parabolic  equation 
model  bottom  prediction  may  be  scaled  to  it  approximately  by  varying  a to 
compensate  for  the  choice  of  g in  the  following  manner. 

Table  1 covers  a typical  range  of  sub-bottom  gradient  values  and  experi- 
mental bottom  loss  values.  Say,  for  example,  that  the  measured  bottom  loss 
at  a grazing  angle  of  20°  is  5 dB.  If  we  suspect  the  bottom  to  be  soft, 
as  would  be  the  case  for  thick  sediment,  then  choosing  g = 1 s-  we  would 
have  to  adopt  a value  near  a = ,005  dB/m  to  achieve  the  appropriate  experi- 
mental bottom  loss.  To  realise  the  same  bottom  loss  at  20°  for  a hard 
bottom,  the  values  g = 10,  a = .05  would  be  adopted. 


TABLE  1.  RELATIONSHIP  BETWEEN  SUB-BOTTOM  GRADIENT, 

VOLUME  ABSORPTION  COEFFICIENT  AND  BOTTOM  LOSS 


Absorption 
coefficient 
a (dB/m) 

Sub-bottom  gradient  g(  sec  1 ) 
(bottom  hardness  parameter) 

Soft 

1 

Medium 

5 

Hard 

10 

Bottom  loss 
at  20° 
(dB/bounce) 

2 

.002 

.01 

.02 

5 

.005 

.025 

.05 

10 

.01 

.05 

.1 

15 

.015 

.075 

.15 

4.3  Choice  of  range  step  size 

As  a result  of  the  approximate  nature  of  the  parabolic  equation  technique 
a number  of  errors  will  be  introduced  into  the  field  solution  as  it  progresses 
in  range  (see  references  5 and  6).  The  maximum  permissible  range  step 
between  two  successive  marching  solutions  is  determined  by  an  error  term 

3n  2 

E range  step  oc  (kA  r 


The  higher  the  vertical  sound  speed  gradient  or  the  higher  the  frequency, 
the  smaller  the  range  step  must  be  to  keep  this  error  suitably  small. 
Consideration  of  the  square  law  dependence  of  this  error  shows  that  it  will 
accumulate  far  more  slowly  over  the  same  total  distance  if  taken  as  the  sum 
of  a large  number  of  small  steps  rather  than  a few  large  ones. 

Choice  of  the  maximum  permissible  step  size  comes  from  experience.  By 
running  the  program  with  a large  step  size  and  again  with  a reduced  step 
size,  the  step  size  will  be  found  to  be  adequate  when  further  reduction 
causes  no  further  change  to  the  predictions.  At  present  the  program  is 
limited  to  1100  steps,  although  this  may  be  increased  by  altering  array 
dimensions  on  any  large  computer.  Program  running  time  is  independent  of 
step  length  as  such,  but  is  proportional  to  the  total  number  of  marching 
solution  steps  involved  in  the  prediction.  To  minimize  computing  costs  the 
range  step  length  should  be  made  no  smaller  than  necessary. 
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Table  2 gives  examples  of  some  step  lengths  found  to  be  satisfactory,  and 
■ some  typical  running  times.  These  running  times  are  generally  well  under 

| half  those  expected  for  similar  runs  using  the  NUSC  program(ref .6) , and  are 

attributable  to  the  use  of  the  simple  starting  field  algorithm,  rather  than 
the  elaborate  normal  mode  procedure  employed  in  the  NUSC  model. 


TABLE  2.  SAMPLE  MODEL  PARAMETERS  FOR  A RANGE  OF  CASES 
SHOWING  COMPARATIVE  RUNNING  TIMES 


Water 

depth 

(m) 

Freq 

(Hz) 

Maximum 
profile 
gradient 
( sec-1 ) 

FFT 

exponent 

ar 

(m) 

Total 
No.  of 
range 
steps 

Running 

time 

(mins) 

4700 

100 

.1 

12 

250 

500 

5.8 

2000 

100 

4 

9 

50 

800 

1 .3 

100 

200 

100 

8 

5 

4000 

4.4 

5.  SOUND  FIELD  INTENSITY  DIAGRAMS 

In  addition  to  the  tables  and  plots  of  propagation  loss  for  one  source  depth 
and  up  to  five  receiver  depths  generated  by  the  program  at  the  conclusion  of  the 
run,  this  version  of  the  program  also  produces  a sound  field  intensity  diagram 
on  the  line  printer  during  the  course  of  the  run. 

The  water  column  and  bottom  is  represented  by  128  line  printer  character 
positions  across  the  page  with  the  surface  at  the  right  hand  end  of  the  carriage. 
In  practice  the  first  32  points  depicting  the  deeper  half  of  the  bottom  are 
suppressed  because  they  usually  contain  little  information.  Since  the  128 
points  must  represent  an  array  of  up  to  2048  points,  each  represents  an  average 
intensity  over  the  appropriate  members  of  the  full  array. 

To  avoid  the  generation  of  diagrams  over  10  pages  long  when  the  full  1100 

range  steps  are  used,  the  user  can  specify  the  approximate  maximum  number  of 

pages  which  the  diagram  should  occupy.  Where  the  total  number  of  lines  expected 

to  be  generated  would  significantly  exceed  this  space,  then  the  program  will 
print  out  every  second,  or  every  fourth  line  or  so  on,  to  condense  the  full 
diagram  to  a more  acceptable  size. 

Different  intensities  are  produced  by  employing  different  printer  characters, 
with  up  to  three  levels  of  overprinting  used.  The  user  can  choose  the  lowest 
value  of  propagation  loss,  and  the  number  of  decibels  per  intensity  increment 
to  give  the  most  useful  diagram.  In  addition  to  white,  ten  shades  of  grey  are 
available.  Increments  of  2 or  3 dB  will  thus  cover  a full  20  or  30  dB  range, 
which  is  generally  adequate  to  define  all  the  features  of  interest  in  the  sound 
field. 

Since  the  field  at  all  depth  points  is  calculated  at  each  step  of  the  marching 
solution,  very  little  extra  computing  is  necessary  to  produce  the  line  printer 
intensity  diagrams.  When  such  diagrams  are  used  in  conjunction  with  the  usual 
single  receiver  depth  plots,  they  often  contribute  materially  to  a better  under- 
standing of  the  sound  field. 


6.  A RANGE  DEPENDENT  TEST  CASE 

To  demonstrate  the  capability  of  the  model  in  predicting  propagation  in  a range 
dependent  situation,  we  consider  the  simple  case  of  low  frequency  propagation  in 
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ar. idealized  shallow  water  duct.  Sound  speed  is  constant  throughout  the  duct, 
and  the  duct  bottom  is  represented  by  a linear  sound  speed  gradient  of  100  s_1  , 
to  provide  a close  approximation  to  a totally  reflecting  hard  bottom.  For  the 
simple  case  of  a constant  in-water  sound  speed,  the  approximations  inherent  in 
the  parabolic  equation  technique  are  minimized  and  the  behaviour  of  the  sound 
field  within  the  duct  is  readily  interpre table  Jay  consideration  in  terms  of  mode 
structure  or  rays. 

Three  cases  are  shown  in  figure  8,  of  the  sound  field  in  a 100  m duct  for  a 
200  Hz  source  at  50  m depth.  In  the  flat  bottom  case  (figure  8(b)),  the  struc- 
ture which  at  first  glance  appears  to  be  bundles  of  rays  coming  from  the  surface 
and  bottom  near  the  source  is  not  rays  at  all,  but  an  almost  symmetric  modal 
interference  pattern.  A maximum  intensity  spot  occurs  at  regular  intervals 
1 .44  km  apart  at  source  depth  with  strong  minima  near  the  surface  and  bottom  at 
the  same  ranges.  At  ranges  halfway  between  these  points,  nulls  may  be  observed 
at  source  depth. 

In  the  second  case  (figure  8(a)),  a linear  sloping  bottom  from  100  m at  the 
source  to  1 50  m at  5 km  range  produces  n similar  pattern,  but  stretched  in  range 
towards  the  deep  end.  The  effect  can  be  even  more  clearly  seen  in  figure  8(c), 
where  the  bottom  slopes  up  from  100  m to  50  m over  the  same  range,  producing  a 
sound  field  pattern  increasingly  compressed  in  range  towards  the  shallow  end. 

This  is  made  clearly  apparent  by  comparing  the  strong  similarities  in  form  of  the 
interference  structure  of  figures  8(b)  and  8(c),  and  noting  the  relative  positions 
of  successive  nulls. 

In  terms  of  ray  theory,  every  reflection  from  the  bottom  in  the  case  of 
figure  8(c)  will  be  at  an  angle  1.14°  greater  then  the  incident  grazing  angle,  so 
that  where  the  ray  equivalent  for  the  lowest  order  mode  is  initially  about  4 , it 
can  be  seen  that  the  effect  of  several  bounces  will  be  quite  significant. 
Calculation  of  the  variation  of  ray  equivalent  cycle  distances  (the  ray  equivalent 
of  a particular  mode  is  the  ray  which  goes  through  one  full  cycle  in  the  same 
distance  as  the  mode  also  goes  through  one  full  cycle)  for  the  main  source  depth 
reinforcements  is  consistent  with  the  range  compression  of  the  pattern  due  to 
the  influence  of  the  sloping  bottom. 

Thus  despite  the  model  being  totally  digital  in  form,  so  that  the  bottom  is  in 
fact  represented  as  a number  of  small  steps  rather  than  a true  slope,  the  model 
can  nevertheless  effectively  simulate  the  reflection  process  at  a sloping  bottom. 
The  reason  is  that  in  the  example  given,  the  bottom  model  is  well  within  the 
diffraction  limits  necessary  to  represent  the  case  as  reflection. 

Figure  9 shows  the  sound  field  for  a case  similar  to  those  presented  in  figure 
8,  but  this  time  with  the  bottom  depth  varying  in  a moderately  irregular  manner, 
rather  than  the  uniform  fashion  of  figure  8.  The  flat  bottom  case  of  figure  8(b) 
is  included  ngain  for  ease  of  comparison.  We  see  that  whereas  in  the  cases  of 
figure  8,  in  which  the  modal  interference  pattern  was  simply  compressed  or  extend- 
ed because  of  the  special  conditions  of  constant  bottom  slope  and  constant  in- 
water sound  speed,  in  figure  9 the  sound  field  intensity  pattern  is  changed 
completely  by  the  irregular  bottom  variation.  Figure  10  represents  curves  of 
propagation  loss  versus  range  for  three  different  receiver  depths  for  the  case 
of  figure  9.  In  addition  to  tables  of  propagation  loss  values,  up  to  five  such 
curves,  one  for  each  receiver  depth,  are  generated  routinely  at  the  conclusion  of 
each  run. 


7.  CONCLUSION 

The  parabolic  equation  model  described  in  this  report  is  an  extremely  versatile 
aid  to  the  prediction  of  c.w.  propagation  loss  prediction  at  low  frequencies  in 
the  ocean.  One  of  the  features  of  the  model  which  perhaps  sets  it  apart  from 
most  other  models  is  the  way  in  which  use  of  it  engenders  a feeling  for  the  inter- 
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action  of  the  sound  field  and  its  environment.  It  is  true  that  a great  deal  of 
insight  into  the  behaviour  of  a sound  field  can  be  gleaned  from  ray  theory 
considerations,  but  at  low  frequencies  attempts  to  visualise  the  field  in  terms 
of  wave  behaviour  are  less  likely  to  be  misleading. 

Since  this  model  could  in  principle  accept  a new  sound  speed  profile  and  a 
new  bottom  depth  at  every  range  step,  the  capacity  of  the  model  to  be  influenced 
by  the  input  data  is  orders  of  magnitude  ahead  of  our  ability  to  provide  such  a 
data  density. 

Such  a situation  highlights  the  need  for  clear  thinking  on  the  part  of  potent- 
ial users  of  propagation  loss  predictions.  When  it  is  reasonable  to  assume  that 
the  full  amount  of  environmental  information  will  never  be  available  and  would 
rarely  be  used  even  if  it  were,  it  is  proper  to  question  the  accuracy  it  is 
reasonable  to  expect  out  of  such  a prediction  model.  Further  it  is  proper  to 
question  the  need  for  developing  a model  free  of  the  limiting  approximations  of 
the  parabolic  equation  method,  when  the  predictions  of  such  a model  may  be  no 
more  accurate  then  those  of  the  parabolic  equation  model  due  to  the  limited 
accuracy  of  the  environmental  data. 

There  are  nevertheless  areas  where  the  parabolic  equation  technique  runs  into 
difficulties,  and  a more  rigorous  approach  seems  the  only  way  on.  The  handling 
of  sloping  bottom  c-ses  where  there  is  a large  and  systematic  change  in  bottom 
depth  in  the  presence  of  oceanic  sound  speed  profiles  is  difficult  for  the  para- 
bolic equation  model,  and  the  modelling  of  an  oceanic  front  also  present  problems. 
Any  model  which  could  hope  to  address  such  problems  with  any  semblance  of  rigour, 
would  presumably  have  to  dispense  with  the  use  of  the  FFT,  and  could  thus  be 
expected  to  generate  prohibitive  amounts  of  computing.  An  investigation  of  this 
problem  being  conducted  by  the  author  at  the  present  time  serves  to  maintain  an 
awareness  of  the  wonderful  power  and  simplicity  of  the  parabolic  equation  tech- 
nique . 

It  has  been  jestingly  said  that  regardless  of  the  prior  work  carried  out  in 
designing  and  testing  any  model,  the  first  time  it  is  ever  used  to  assist. in 
someone  else's  problem,  the  case  presented  turns  out  to  be  either  one  which  finds 
a residual  error  in  the  program,  or  else  one  quite  beyond  the  scope  of  the  model. 
The  presently  described  model  is  no  exception,  since  there  are  limits  to  its 
applicability.  It  is  hoped,  however,  that  this  report  will  make  the  limits 
quite  clear  to  any  potential  user,  and  that  the  model  is  sufficiently  simple  to 
operate  that  with  a little  experience  any  new  user  may  generate  results  which 
instil  confidence. 
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APPENDIX  I 

INPUT  DATAIOR  THE  PARABOLIC  EQUATION  MODEL 

1.1  All  data  (unless  alphanumeric)  is  read  in  under  formats  110  and  FI 0.0, 

which  means  that  all  whole  numbers  are  right  justified  to  character  position 
10,  20,  30  etc.,  and  any  numbers  with  decimal  point  included  may  occur  any- 
where in  the  appropriate  floating  point  field. 


DATA  DECK  STRUCTURE 


Block 

Card 

Data 

1 

1 

TLPLOT,  MFLAG 

2 

1 

RMAX,  DR 

3 

1 

RSCALE,  TLMAX 

4 

1 

SIGMIN,  SIGSCL,  PAGES 

5 

1 

BETA 

6 

1 

SD,ANG 

7 

1 

nr,(rd(ii),  II  = 1 ,NR) 

TITLE  (Alphanumeric,  60  characters) 

8 

1 

9 

1 

H,C,GRAD 

10 

1 

F 

11 

1 

MFFT,NSEC 

12 

1 

PRANGE ,BB0T 

onwards 

2 

ND 

3 

ZL(1 ),cl(i ) 

(NSEC 

Blocks) 

zl(nd),cl(nd) 

i 


1.2  Explanation  of  input  data 


TLPLOT 

MFLAG 

RMAX 

DR 

RSCALE 

TLMAX 

SIGMIN 

SIGSCL 

PAGES 


BETA 

SD 

ANG 

NR 

RD(II) 


plot  flag  set  to  1 for  plot;  0 for  no  plot 
1 for  metric  units  plot;  0 for  English  units 
maximum  range  for  transmission  loss  calculation  (m) 
range  increment  (m) 

desired  plot  range  scale  in  kyds/inch  (English)  or 
km/inch  (metric) 

maximum  value  of  transmission  loss  axis  (scale  fixed 
at  10  dB  per  inch) 

minimum  value  of  propagation  loss  for  grey  scale  of 
line  printer  sound  field  intensity  diagrams 
number  of  decibels  per  grey  level  step  for  sound  field 
intensity  diagrams 

length  control  for  intensity  diagrams.  Number  of 
printer  pages  will  lie  between  one  and  two  times  this 
value,  up  to  a maximum  correponding  to  one  line  per 
range  step 

sub-bottom  attenuation  (dB/metre) 
source  depth  (m) 

limiting  grazing  angle  of  interest  at  the  source 
number  of  receiver  depths  (maximum  5) 
depth  of  (il)th  receiver 
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TITLE 

H 


c 


GRAD 

F 

MFFT 

NSEC 

PRANGE 

HBOT 

ND 

ZL(l  ),  CL(1 ) 
etc 


header  for  transmission  loss  plots;  centre  in  first 
60  characters 

arbitrary  water  depth,  usually  chosen  to  be  that  at  the 
source.  This  value  scales  the  whole  propagating  field 
to  be  4 x H (m) 

uniform  speed  of  sound  in  bottom  below  sub-bottom 
gradient  layer  (m/sec) 

sub-bottom  sound  speed  gradient  (s-^) 
frequency  (Hz) 

2**MFFT  is  size  of  the  Fourier  transform 

must  be  greater  than  or  equal  to  2.  Total  number  of 

sound  speed  profiles  to  be  read  in. 

distance  from  source  in  metres  of  associated  sound  speed 
profile 

bottom  depth  (m)  applicable  to  PRANGE 

number  of  points  on  sound  speed  profile  for  individual 

range  block 

depth,  sound  speed  pairs  in  metres,  metres/sec.  Depth 
is  zero  at  sea  surface,  and  direction  of  positive  increase 
is  downwards. 


BOTTOM  LOSS  (dB) 


Figure  2.  Bottom  loss  versus  grazing  angle  as  a function  of 
processing  bandwidth  (center  frequency  = 110  Hz) 
(from  ref.10) 
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(a)  Sound  speed  in  water  1500  m.s-1,  bottom  sound  speed  1550  m.s" 
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(b)  Sound  speed  in  water  1500  m.s  , bottom  sound  speed  1700  m.s 


Figure  3«  Rayleigh  reflection  loss  curves  for  a number  of  bottom 
densities  and  two  bottom  sound  speeds 


BOTTOM  LOSS  (dB) 


Figure  4.  Theoretical  bottom  loss  for  the  Tyrrhenian  Sea  (from  ref. 9) 


20 
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3ce  of  the  ray-equi 
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Figure  8.  Sound  field  in  hard  bottom  duct  for  200  Hz  source  at  50  m depth 


PROPAGATION  LOSS  (dB/lm) 


) 


Figure  10.  Propagation  loss  versus  range  curves  for  3 

receiver  depths  (ZR)  for  the  case  depicted  by 
figure  9 
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